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WL-TO-61*-10l* 


ABSTRACT 


A  nuaber  of  commonly  used  finite-difference  analogs  to  partial  derivatives 
in  tvo  space  distensions  are  investigated,  and  a  few  variations  are  proposed. 
The  accuracy  of  these  analogs  is  assessed  by  obtaining  numerical  results  for 
deformations  for  which  the  analytical  gradients  can  be  evaluated.  None  of 
the  analogs  appeared  superior  for  those  deformations  which  were  investigated, 
and  it  appears  that  a  choice  may  be  made  on  the  basis  of  computational  con¬ 
venience. 
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SECTION  I 


INTRODUCTION 


One  of  the  principal  difficulties  in  constructing 
finite  difference  analogs  to  the  Lagrangian  equations  of 
compressible  fluid  motion  in  two  space  dimensions  is  the 
choice  of  finite  difference  representation  of  the  spatial 
derivatives  in  the  momentum  equation.  A  number  of  different 
expressions  have  appeared  in  the  literature,  but  there  seems 
to  be  little  information  available  regarding  their  relative 
merits. 

Finite  difference  representations  of  spatial  gra¬ 
dients  should  introduce  minimal  truncation  errors.  However, 
some  other  properties  are  often  considered  to  be  more  im¬ 
portant.  As  a  Lagrangian  mesh  becomes  severely  distorted, 
it  is  found  that  some  finite  difference  expressions  lead  to 
accelerations  in  the  wrong  direction.  This  leads  to  an  in¬ 
stability  and  eventual  stoppage  of  the  computer.  Such  se¬ 
vere  distortions  sometimes  occur  in  relatively  unimportant 
areas  of  the  flow,  and  it  is  often  considered  more  impor¬ 
tant  that  a  finite  difference  expression  lead  to  accelera¬ 
tions  of  the  correct  sign  than  of  the  correct  absolute  mag¬ 
nitude. 

It  has  occasionally  been  argued  that  finite  dif¬ 
ference  expressions  which  show  large  truncation  errors  lead 
to  smoother  solutions,  and  are  preferable  on  this  account. 

It  may  be  noted  that  truncation  errors  sen/e  to  increase 
the  acceleration  in  areas  of  high  gradient.  In  this  re¬ 
spect  truncation  operates  in  precisely  the  same  way  as  ar¬ 
tificial  viscosity,  which  is  formulated  to  provide  an  ad¬ 
ditional  acceleration  in  areas  of  high  gradient.  It  seems 
preferable,  however,  to  choose  the  finite  difference 
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expression  which  has  least  truncation  error,  leaving  smooth¬ 
ing  of  the  solution  to  artificial  viscosity  terms  over  which 
one  can  exercise  some  control. 

Several  common  finite  difference  expressions  are 
discussed.  In  order  to  investigate  truncation  errors  the 
following  procedure  was  followed.  A  quadratic  deformation 
was  postulated,  for  which  the  density  change  could  be  found 
analytically  using  the  principle  of  conservation  of  mass.  On 
assuming  an  equation  of  state,  the  associated  pressure  change 
could  be  found  analytically,  and  the  acceleration  found  by 
applying  the  principle  of  conservation  of  momentum.  A 
Murnaghan  equation  of  state  was  used.  The  corresponding  de¬ 
formation  of  an  initially  square  mesh  was  calculated  for  the 
same  quadratic  deformation,  so  that  the  acceleration  could  be 
computed  using  the  finite  difference  expressions.  The  be¬ 
havior  of  the  finite  difference  expressions  after  crossing 
of  cell  vertices  and  after  cell  inversion  was  also  investigated. 
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SECTION  II 


DIFFERENTIAL  EQUATIONS 


2.1  Tensor  Equations 

Denote  the  spatial  (Eulerian)  coordinate  and 
the  material  (Lagrangian)  coordinate  Xx,  where  we  use  majus¬ 
cules  to  refer  to  the  original  state,  minuscules  to  refer  to 

•  «*■  _ 

the  current  deformed  state.  At  time  t“0  take  x'  ..  X 
Consider  the  deformation 

X  i  .  x;  (  X',  X‘  ,  X‘)  2.1 


with  Jacobian  J  given  by 

t  -  I  *  ;r  1  2.2 

where  (  ),j  is  the  covariant  derivative  with  respect  to  the 

Xx. 


The  tensor  equation  of  mass  conservation 

£•  .  T  -  — 

P  dV 


is  simply 

2.3 


where  f  is  the  density,  pa  the  initial  density  at  time  t=0  and 
di r  and  d V  corresponding  volume  elements  in  the  deformed  and 
undeformed  states.  The  tensor  equation  of  momentum  conservation 
in  a  perfect  fluid  is 

*  -  p  P,i  2.4 

where  a1  is  the  acceleration  vector  , /»  the  pressure  and  ( 
is  the  covariant  derivative  with  respect  to  the  x* .  It  will 
be  useful  to  write  this  equation  in  the  equivalent  form 

<*i  *  -  jfc  P,  I  2.5 


or  since 


CO  fait or 


T 


2.6 


3 


2.  / 


aL  *  '  %  x  \z  ) 

where  we  have  used  equation  2.3. 

We  will  choose  a  quadratic  deformation  (equation  2.1) 
in  the  form 

'  <*-T  X'  -  bi  2.8 

where  the  coefficients  c  *  ^-x  and  A*  are  arbitrarily 
chosen  constants.  It  can  be  seen  that  the  b*  introduce  a  uni¬ 
form  translation  while  the  introduce  a  homogeneous  deforma¬ 
tion  and  rotation.  These  terms  do  not  lead  to  accelerations 

i 

and  are  retained  for  convenience  in  controlling  the  size  and 
position  of  the  deformed  meshes.  Since  we  have  not  restricted 
ourselves  to  small  deformations  the  and  cannot  be  in¬ 

terpreted  easily  as  rotations  and  deformations,  but  it  may  be 
noted  qualitatively  that  the  diagonal  terms  introduce  stretches 
in  the  corresponding  coordinate  directions,  while  the  symmetric 
parts  of  the  off-diagonal  terms  introduce  shears,  the  anti¬ 
symmetric  parts  rotations.  The  formulation  of  the  deformation 
in  equation  2.8  is  therefore  quite  flexible  in  producing  almost 
any  desired  deformed  mesh  shape  by  appropriate  choice  of  con¬ 
stants. 

Once  the  deformation  (equation  2.8)  has  been  speci¬ 
fied,  the  density  distribution  may  be  found  from  equation  2.3. 
It  is  necessary  to  postulate  a  relation  between  p  and  p  before 
finding  the  acceleration.  The  particular  form  of  this  relation 
does  not  affect  the  conclusions  materially.  An  "equation  of 
state"  in  the  Murnaghan  form  has  been  chosen  for  its  simplicity 

p  •  k  {(£)*  ■ 1 } 

or  using  equation  2.3 

p  -  k  (  J"r  “  I  )  2.9 
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where  k  and  }f  are  material  constants.  We  are  principally  con¬ 
cerned  with  relatively  incompressible  materials,  and  constants 
appropriate  for  aluminum  have  been  used  in  the  calculations  de¬ 
scribed  later  (i.e.  k*  /  1  *  10"  dyn  cm^  ,  3r*  4  ,  f,  s  2-7P  gm  cm  ^)  . 


2.2  Physical  Equations 

Only  two-dimensional  rectangular  Cartesian  coordinates 
(  x,  x  )  and  two-dimensional  axially  symmetric  coordinates 
(  x,  x  where  we  have  written  x  for  the  radius)  are  considered. 
The  continuity  equation  (equation  2.3)  is  usually  taken  in  its 
second  form 

p  djr 

Po  S  <*v  2.10 

while  the  momentum  equation  (equation  2.4)  takes  the  form 

i 

a*  =  "  jo 


i  b? 
«X  -  -  J  fZ 


2.11 


Defining  «*/  for  the  rectangular  case,  and  « »  2  for 
the  cylindrical  case,  the  Jacobian  for  both  cases  may  be  written 
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d  X 

bZ 


bx 

bl 


2.12 


Thus  the  alternate  form  of  the  momentum  equation 
(eouation  2.7)  becomes 


dp  i*  )  /  *  )“ 

di  9x  J  '  x  J 


2.13 


a 


±  (  i£  lr  _ 
fi,  (  dX  dZ 


l  f  ip  ix  d£  dx  )  <  *.  I-'1 
4*  *  '/>.  1  iX  "  dX  dZ  /(  X  / 

2.3  Analytical  Gradients 

Specializing  to  two-dimensional  rectangular  Cartesion 
and  cylindrical  polar  coordinates,  the  deformation  transforma¬ 
tion  (equation  2.8)  retains  only  the  following  terms 

x  *  cn  X 1  *  ctx  X  £  ♦  x 

♦  d„  X  *  d„  Z  *  b, 

2.14 

2  .  X  ♦  c.j  X2  ♦  c,j  2 

+  <*„X  *  d„  Z  *  4j 


Thus  the  Jacobian  of  the  transformation  can  be  writ¬ 
ten  (from  equation  2.12) 

7  *  (A  D  -  BC)  £""  2.15 

where 


a  2  f  y 
/*  ■  *  " 


* 


a  M 


2  c,9Z  +  c,x  X  +  A, % 
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C  --  fa  *  2  £JI  x  *  *  *  d*' 

d  “  M  * 2 

E  x  C„  X  +  C,\  2  +  d„  +  (c,9  2l  +  d,3  2  +  b,)/X 

Then  density  and  pressure  are  given  bv  equations  2.3 
and  2.9  respectively.  In  order  to  find  the  accelerations,  it  is 
necessary  to  use  the  second  form  of  the  momentum  equations, 
which  may  be  written 


Note  that  differentiation  of  equation  2.9  leads  to 


if 

3X 


-  *  k 


J 


aj 

ax 


aP 

32 


-  y  k  J 


-(*+') 


37 

32 


so  that  equations  2.16  can  be  conveniently  written 


2.17 


a 


K 


a 


z 


m  i  b.  E  *  ' 
%  j  Tr> 

f  k 
P. 


2.18 


7 


where  the  derivatives  of  J  are  easily  found  by  differentiating 
equation  2.15. 

Once  the  values  of  the  12  coefficients  are  specified 
in  the  deformation  transformation  (equation  2.14)  it  is  a 
straightforward  task  to  evaluate  the  density,  pressure  and  ac¬ 
celerations  at  any  particular  point  specified  by  its  initial  co¬ 
ordinates  (  X,  2  ) . 

Finally  the  magnitude  of  the  acceleration  vector  is 

*  *  / 2.19 

and  the  angle  between  the  acceleration  vector  and  the  x  axis 
is 
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SECTION  III 


FINITE  DIFFERENCE  EXPRESSIONS 

A  number  of  finite-difference  analogs  to  the  equation 
of  mass  and  momentum  conservation  have  been  evolved.  An  attempt 
is  made  to  collect  brief  derivations  of  some  of  these  here. 

Quantities  are  considered  only  at  a  finite  number  of 
locations  in  space,  initially  distances  4X  and  42  apart.  The 
initial  X  coordinate  after  the  I  *  increment  AX  is  denoted  Xx  , 
and  the  initial  2  coordinate  after  the  J *  increment  ^Fis  de¬ 
noted  .  In  effect  the  material  is  covered  by  a  finite  co¬ 
ordinate  grid  which  deforms  with  the  material  (Fig.  1).  Co¬ 
ordinates  x  and  z  at  time  t  for  the  point 
denoted  xZ  J  ,  aI  7. 

While  positions,  velocities  and  accelerations  are  con¬ 
sidered  only  at  the  vertices  of  the  finite  difference  grid,  den¬ 
sities  and  pressures  are  considered  averaged  over  the  meshes, 
and  are  denoted  7+„x  ,  etc. 

In  developing  the  finite  difference  equation,  it  is 
considerably  more  convenient  to  use  the  notation  of  Figure  1, 
translating  the  equation  into  indicial  notation  before  program¬ 
ming  for  the  computer. 

3.1  Mass  Equation 

The  second  form  of  the  mass  equation 

f  ‘  it  3a 

can  be  set  into  finite  difference  form  by  considering  the  de¬ 
formed  mesh  to  consist  of  quadrilaterals.  The  current  area  of 
mesh  1  (Figure  1)  is  then 

A,m  i  {(x*~  *o)(zo~2*)  X*»  -**) J  3,2 


Ai  ’  %  are 
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while  the  original  area  X, 0  at  time  t»0  is  found  by  sub¬ 
stituting  X  and  2  for  x  and  z  in  equation  3.2. 

Both  the  rectangular  Cartesian  and  cylindrical  polar 
cases  may  be  written  down  simultaneously  by  defining  *  / 
for  the  rectangular,  a  •  2  for  the  cylindrical  cases  as  be¬ 
fore.  The  mass  equation  becomes 

F'  '  /».  TzTr™  3,3 

where  x,  is  the  radius  of  the  centroid  of  the  area  /I,  ,  and  m, 
is  a  mesh  constant  defined  as 

-  m  t  r,  1*1 


m,  •  A,'  (  ) 


The  radius  of  the  centroid  appears  only  In  the  cylindrical  case. 
A  reasonable  approximation  for  moderate  distortions  is 


i  ( x«  ♦  *«  ♦  *»  *  *-) 


'  44 

N 

\ 

i.  \ 
_  \ 


A  better,  but  lengthier,  expres¬ 
sion  may  be  obtained  by  dividing 
the  quadrilateral  into  two  tri¬ 
angles,  i.e.,by  taking 


( Si)*"1  +  A*  ( *»)" 


where 


i  {(**-**/*• -*J 

a  •  i  { (*•  - *4  **  -*•)-(*>-  *4**  -*•)  J 


•  %  (  **  ♦  **  ♦  *,  ) 


xc-  %  (  x9  ♦  x# 


♦  X* 
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3.  7 


and  where  m ,  is  given  by 

m.  *  A  (  d.*  f  V  j*  '  ♦  (*m‘)  "  '  J 

For  the  rectangular  case,  *•/,  and  these  equations 
are  identical  to  the  previous  ones,  (equations  3.3  and  3.4). 

It  may  be  noted  parenthetically,  that  the  same  results 
are  achieved  by  starting  with  the  alternate  form  of  the  mass 
equation  (equation  2.3)  as  may  be  expected.  Terms  in  the  Jaco¬ 
bian  (equation  2.12)  may  be  written  in  finite  difference  form  as 


» 

■  » 

/  **  +*H  _  *o  +  ) 

ax 

\  X  X  ' 

a* 

-  —  / 

f  *»  +  **  ) 

3.8 

a* 

42  \ 

X  X  j 

etc . 

When  these  are  substituted  into  equation  2.12,  and  terms  are  ex¬ 
panded  and  simplified,  we  obtain 

J-  >  '  {(Xx-X'KZ‘-Z»)  -  (**-*»Xx’-**J}  3.9 

When  this  is  inserted  into  equation  2.3,  we  immediately  arrive 
again  at  equations  3.3  and  3.4. 

3.2  Momentum  Equation 

A  variety  of  difference  analogs  have  been  given  for 
the  pressure  gradient  terms  in  the  momentum  equations  (equations 
2.11).  Some  of  these  are  described  in  the  following. 


3.2.1  Taylor's  Expansion 

The  pressure  is  known  at  discrete 
points  1,  2,  3,  4  surrounding 
point  0  at  which  the  pressure 
gradients  are  to  be  evaluated. 
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One  method  is  to  apply  Taylor's  expansion  between  points  0  and 
1,  0  and  2, etc.  to  obtain  four  equations  of  the  fomi. 


•  p.  *  (*■  -*•)  H  *  (x-  -*•)  ii 


•  I  „  ,  yp  ,  / 

+  x(*r”*)  *'-**) jTTx  *  *  (*•  JI*  r  "" 


3.10 


where  terms  of  third  and  higher  order  in  (  *,-*•)  and  (x,-*0) 
have  been  omitted.  Kolsky ‘proposed  solving  this  overdetermined 
system  of  equations  for  and  ££  by  first  solving  for  (P, -^) 
and  (A.-/U).  obtaining  the  two  equations 


*  H (*''**) 


*  ik  (*'•**)*. 


*  *  »t£x [(*•  s*'>  *  (■  *' 


{PX~P*)  ’  $X  (**’*+)  *  Jx  (z'~z*) 

*  *  Yzx  (*L~X+) 

*  *  2 firfz  ~  ( Xx  ~**)  *  " 

where  SMn  *  f  (  *»  ♦  **  )  -  *•  ,  etc.  are  measures  of  the  assym- 
metry  of  the  mesh.  Omitted  terms  are  of  second  and  higher  order 
in  the  mesh  size  (  ■*.  ** *> ) »  etc .  Providing  the  mesh  size  is 
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small,  these  higher  order  terms  are  negligible  compared  to  the 
terms  retained  above.  The  above  equations  may  now  be  solved 
for  the  gradients 


■  j»  /*  *  ■  rxj  -  (f>,  ~2*)}  *■  /f« 

3.11 


where 


*  i  **•**)  - 1  *>-*+)} 


It  is  ceen  that  ft  represents  the  area  of  the  quadrilateral 
12  3  4. 

The  remainder  terms  ft*  and  ftx  involve  products  of 
the  mesh  size  and  the  £ '*  (e.g.,  (* ,  -  X  )  (f  x  tJ  ,  etc.).  If  the 
mesh  is  nearly  symmetric,  these  terms  are  negligible  compared 
to  the  terms  which  have  been  retained.  It  is  clear  that  the 
error  terms  vanish  for  an  undistorted  mesh,  but  become  progres¬ 
sively  larger  as  the  mesh  distorts. 

Equations  3.9  require  coordinates  of  the  centroids  of 
quadrilateral  1,  2,  3,  4.  Kolsky  simply  used  the  approximate 
equation  3.5,  obtaining  the  final  equations 

m  tUif)  [Ia-hX*'-**  *  *'  'Zh) 

~(p.-p,X  *<  -  ^  *•  •*«  -**J) 


a*‘  tfrri 


■f(p>-p,)(x<  -*»  +  -*») 

-  (p,-  P,x  *<  -  x»  +  *a  *  *<■  -  *r)} 


3.12 


n 


There  are  two  possibilities  for  representing  the  denominator 
(*/>  )  in  equation  3.12  in  a  simple  form.Kolsky  used1 

(»fi)  *  ?  (  *'P,  *  *  *1*)  3*13 

A  second  expression  follows  from  equation  3.3  for  the  cylindri¬ 
cal  case 

(/ip)  *  i  (  ",  +  "x  ♦  "i  ♦  m+  )/(*•)"'*  3.14 

where  is  used  as  an  approximation  to  the  centroid  of  the 
quadrilateral  1  2  3  4.  Equation  3.7  is  used  for  the  m's .  For 
certain  serious  distortions  this  may  lead  to  considerable  error. 
For  the  rectangular  case,  equation  3.14  is  equivalent  to  equa¬ 
tion  3.13. 


It  may  be  noted  parenthetically  that  the  same  results 
are  achieved  by  starting  with  the  alternate  form  of  the  momentum 
equation  (equation  2.13),  as  may  be  expected.  The  Lagrangian 
pressure  gradients  may  be  represented  by 

_L  /  A  ♦  A  _  P»  ,  ,e 

AX  (  X 


3X 


IP 


AX 

I 

az 


Px  +Pj 


) 


Inserting  these  together  with  equations  3.8  into  equations  2.13 
and  simplifying, 


‘■•'4r£rz7 

■■  •  (<*■'** "•>) 


3.16 


With  the  aid  of  equations  2.3  and  3.9,  (the  latter  written  for 
quadrilateral  1234),  it  is  seen  that  the  above  equations  3.16 
are  exactly  equivalent  to  the  previous  result,  equations  3.11. 
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3.2.2  "Midpoint"  Method 

Amurud  and  Orr^  noted  that  Kol- 
sky's  scheme  led  to  reversal  of 
signs  of  the  accelerations  when 
the  mesh  became  sufficiently  dis¬ 
torted  so  that  quadrilateral 
1234  became  inverted.  They 
proposed  applying  Taylor's  theo¬ 
rem  between  points  0-5,  0-6,  0-7 
0-8  as  shown,  in  conjunction 
with  a  test  and  approximate  cor¬ 
rection  procedure  to  prevent  in¬ 
version  of  quadrilateral  5678. 

Substituting  subscripts  5,  6,  7,  8  for  subscripts  1, 

2,  3,  4  in  equations  3.11  and  writing 

ps  ■  i  (  *  P\)  etc. 

m  i  (  **  +  *•)  etc* 

the  following  equations  result 

-*«  ♦  -  *•)} 


A 


•  / 
C/o 

-  J 

7  4 

c 


3.17 


The  same  expressions  may  be  used  for  (flf)  as  before  (equations 
3.13  and  3.14). 


3.2.3  Green’s  Transformation 


A  different  approach  to  the 
problem  of  finding  pressure 
gradients  follows  from  Green's 
Transformat ionf'*  In  two  dimen¬ 
sions 


where  nL  is  the  outward  normal  vector  to  the  circuit  f  enclos¬ 
ing  the  area  ft  ,  and  the  comma  denotes  covariant  differentia¬ 
tion.  In  component  form,  we  are  led  to  the  two  equations 


/ 


Pti  Ack 


3.18 


i**'  *  f%  **  -  (kL  * 

-jpix  -  dfi  -  (kL„  1 9 


3.19 


Choosing  the  circuit  A  B  C  D,  and  considering  the 
pressure  on  side  A  B  to  be  given  by  its  average  value  px  etc.. 


{*  t**'*’-1  **■(*••**!* 

3.20 


where  ft  is  one  half  the  area  of  the  quadrilateral  A  B  C  D.  A 
reasonable  approximation  to  the  denominator  is  given  by  equation 
3.13  or  equation  3.14.  For  irregular  distortions,  a  better  but 
much  lengthier  expression  would  be 

ftp  m  i  (  ft/p,  ♦  ft*'p*  +  ft*  p>  +  ft+p+  )  3,21 
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where  P, '  is  the  area  of  triangle  A  0  D  etc. 

Variations  of  equations  3.20  suggests  themselves.  In¬ 
stead  of  representing  (  ftp  )  by  equation  3.21,  we  write 


*  -±)  ,  *  (*•  •**)  * 
1  A  Pi  2  *x  px 


.  *(*,-*J  l 

2*iPi  J 

3.22 


f  P>(**:Xq) 
<  2*,P, 


P*(X»  -**)  +  P»  (*€  -*b)  +  P*  (  X,  -Kc)  J 

2  P*.  A.  2  Pi  Pi  2  *lP+  J 


where  A,  etc.  are  the  areas  of  quadrilaterals  A  0  D  H  etc.  In 
these  equations  the  contribution  of  each  mesh  to  the  accelera¬ 
tion  is  independent  of  the  other  meshes.  We  can  also  use  equa¬ 
tion  3.3  to  write 

*>p>  *  (irt®'1  3-23 

etc.,  where  we  use  equation  3.5  for  x,  ,  in  each  of  the  above 
terms.  A  further  variation  is  to  write 


A  P,  -  1  A>,  3.24 

etc.,  where  P/  is  the  area  of  triangle  A  0  D,  in  each  of  the 
above  terms. 


It  may  be  noted  parenthetically  that  if  circuit  123 
4  is  chosen,  the  results  obtained  previously  by  means  of  Tay¬ 
lor's  Theorem  are  obtained.  The  average  pressure  on  side  1  2 
can  be  written  i  [  P >  +  P*  )  .  Thus,  applying  equations  3.19 

to  the  circuit  1  2  3  4, 

«,  •  I (Pf)’, )(**-*■)  *  (Pt  +  P*Xx>  ~*J 

+  +(*  ~p+X*>  -*■>) }  3.25 

’f*T  +(P> 

* (p.+Pi)(x.-*>)  *  (p.+p*X*,-* .)} 
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where  A  is  the  area  of  quadrilateral  1234.  Upon  simplifica¬ 
tion  and  interpolation  of  x  and  x  these  reduce  to  equations 
3.12. 

Green's  Transformation  has  also  been  used  to  obtain 
expressions  of  the  gradients  in  a  triangular  mesh?  Triangular 
meshes  are  generally  found  to  be  anisotropic  and  are  not  con¬ 
sidered  here. 


3.2.4  "Force  Gradients11  Method 

An  expression  for  the  gradients 
in  use  at  L.A.S.L.^may  be  de¬ 
rived  in  several  ways.  For  the 
shaded  zone,  Lagrangian  gradients 
may  be  expressed  as 

/if\  P,-P4  /ii  )  *» 

wJ,  *  *  ITT. 

Writing  similar  expressions  for  the  zone  2  B  3  0  on  the  opposite 
side  of  0,  and  interpolating,  the  X  gradient  at  0  can  be 
written  (from  equation  2.13) 


_i  /  l^_  I"'1  (r*-nXz*  •*•)  i  g  l*'1 

**  *  (  (X.-X4j|f*,-Z.J  1  x/t  +  (Xt-XJZ.-Z,) 

3.26 

II  l‘"  II]-'  ) 

(*♦ -2j[X.-X')[X/7  ) 

with  a  similar  expression  for  Q M.  Two  ways  of  representing  the 
denominators  of  the  terms  above  follow  from  the  mass  equation. 
We  note  that  the  denominator  of  the  first  term  (X,“X4) 
represents  the  undeformed  area  of  the  shaded  rectangle,  denoted 
by  (  )f.  Using  equation  3.4 

f,,(X,-X*X2t>-z°)(x)  ,=  mr  “  *(m>*m*) 
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Proceeding  in  this  way,  the  accelerations  become 


where  we  have  simply  approximated  ae^  *  £  (■*&  • )  etc. 


Alternatively,  using  eouation  3.3  these  equations  can 
be  put  into  the  form 


[p,  -/%X  *D-*o) 

(Pk-PjX*-  -  *■•) 

P,P,  ■>  P+ 

A  *  fi>  Pi 

(r,  -kX*»  -*•) 

(  Ip*- PiX  *•-*<) 

(Px-Pi  X  *•-*•) 

P,*,  *  P*  ** 

P\  Ax  +  ft  As 

[Pt-PiX**  ~*°) 

(  J°4  ~  PsX  ~  *<■) 

P.fl,  *P<  *x 

P*  *  P,  ^ 
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The  two  approximations  made  in  the  denominators  of 
equations  3.27  and  3.28  are  analogous  to  the  approximations 
made  in  enuations  3.14  and  3.13. 

The  above  equations  (3.27  or  3.28)  lead  to  difficul¬ 
ties  for  certain  serious  distortions  in  which  one  or  more  of 
the  segments  OA,  OB,  OC ,  or  00  approach  zero  length.  More  com¬ 
plex  equations,  apparently  derived  on  a  similar  basis,  but  in¬ 
volving  multiple  interpolations  to  overcome  this  drawback,  have 

9 

been  reported . 


3.2.5  Virtual  Work  Method 


1 

■  1 

0 

J 

-t 

Goad^  proposed  a  method  of  obtain¬ 
ing  the  apparent  force,  and  hence 
the  acceleration,  on  a  mass  point 
located  at  0  by  the  principle  of 
virtual  work.  Suppose  that  the 


point  0  is  displaced  by  an  infinitesimal  vector  displacement  Ss k 
Work  is  done  on  the  surrounding  four  cells,  and  the  energy  of 

the  system  is  modified  by  an  increment  $£  given  by 

S  £  *  £  (p  Ss  k  trlh  )i 

where  ir  is  the  cell  volume,  and  the  comma  denotes  covariant  dif¬ 
ferentiation.  Thus  a  force  F.  acts  at  0  given  by 


Thus 


JE  *  F„  Xik 


F*  *  £ 


and  the  acceleration  at  0  is  therefore 

ak  ■  ~  £  ( p  'r,k  ^ 

where  m  is  the  mass  which  is  thought  of  as  being  concentrated 

at  point  0. 

Using  the  second  method  of  finding  the  volume  of  a 
mesh,  (Subsection  3.1)  for  mesh  1  we  have  within  a  factor  (27r>oC ) 


(aj—  ♦  Au(*.r"} 
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where  ,  :*tetc.  are  given  following  equation  3.6.  To  find 
the  derivatives  of  v  ,  we  consider  ,  z0  variable  with  the 
other  cell  vertices  held  fixed,  and  obtain  for  the  x  direction 
for  cell  1 

(f>  ,),  *  »  r>  {  (  *  3  - - ~) 

f  1  Zp  ~2C)( X*  "*•)  ~  (  ~  x»  )J  J 

while  for  the  z  direction  for  cell  1 

i,*..),  ■  i p, { (SLi^r  (*>-**)} 

The  mass  m  associated  with  point  0  may  be  taken  to  be 
(within  a  factor  (2 ir)"'1) 

m  ■  +  (  /  *■  m K  t  tvtj  +  4  J  3.29 

where  m,  etc.  are  given  by  equation  3.7. 

Summing  terms  for  all  four  meshes 

•  A ((*' )"'(*.-*.) * 

♦ J  ^  _  It)  t  ^  -*.)]] 

3.30 

..  •  £  I »  t  *  ft 

. *  ( siiiiff--*.;  -  ft  ( -«•;} 
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For  the  rectangular  case,**/,  and  the  equations  reduce  identi 
cally  to  these  obtained  by  Green's  Transformation  (equations 
3.20). 
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SECTION  IV 


RESULTS  AND  DISCUSSION 


Results  of  a  number  of  calculations  are  given  in  this 
section.  Both  the  analytical  and  finite-difference  equations 
were  programmed  for  computation  on  an  IBM  7094  Computer.  Results 
are  shown  in  the  figures,  where  results  of  the  various  finite 
difference  equations  are  labeled  as  follows. 


1. 

Taylor  Expansion 

Eq . 

2. 

Taylor  Expansion 

Eq. 

3. 

Midpoint  Method 

Eq . 

4. 

Midpoint  Method 

Eq . 

5. 

Greens  Transformation 

Eq. 

6. 

Greens  Transformation 

Eq. 

7. 

Greens  Transformation 

Eq . 

8. 

Greens  Transformation 

Eq . 

9. 

Greens  Transformation 

Eq . 

10. 

Greens  Transformation 

Eq. 

11. 

Force  Gradients 

Eq . 

12. 

Force  Gradients 

Eq. 

13. 

Virtual  Work 

Eq . 

In  each  case  the  accelerations 
nitude  a  and  direction  9  by 


3. 

,12 

with 

Eq  . 

3. 

.13 

3. 

,12 

with 

Eq  . 

3, 

.14 

3. 

,17 

with 

Eq  . 

3, 

,13 

3. 

,17 

with 

Eq  . 

3, 

,14 

3. 

,20 

with 

Eq . 

3. 

,13 

3. 

,20 

with 

Eq . 

3, 

.14 

3. 

,22 

3. 

,22 

with 

Eq . 

3. 

,23 

3. 

,20 

with 

Eq . 

3. 

,21 

3. 

,22 

with 

Eq . 

3. 
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3.28 

3.27 

3.30 

o-*  and  a*  were  converted  to  mag- 


4.1 


0  - 


arc 


tan 


and  percentage  errors  in  magnitude,  and  angular  errors  in  direc¬ 
tion  compared  to  the  true  magnitude  and  direction  were  found. 


It  is  difficult  to  present  results  to  illustrate  the 
growth  of  the  error  as  a  function  of  some  measure  of  deformation. 
Since  the  equation  of  state  (equation  2.9)  was  nonlinear,  the 
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results,  expressed  as  a  percentage  of  true  acceleration,  would 
be  affected  appreciably  by  the  addition  of  a  uniform  (hydrostatic) 
compression.  Moreover, in  the  cylindrically  symmetric  case 
the  results  would  be  heavily  affected  by  a  uniform  translation 
in  the  radial  (x)  direction.  In  order  to  eliminate  these  effects 
to  some  extent  the  following  procedure  was  adopted.  The  quad¬ 
ratic  coefficients  ( C  )  were  adjusted  incrementally  in  some  ar¬ 
bitrary  manner.  The  diagonal  linear  coefficients  (  d,,  , 
were  then  given  values  such  that  the  mesh  with  largest  area  was 
reduced  to  the  area  it  had  before  deformation.  Finally  the  co¬ 
efficient  b ,  was  adjusted  so  that  the  radius  of  the  mesh  vertex 
with  minimum  radius  was  set  equal  to  the  minimum  radius  before 
deformation.  This  procedure  did  not  eliminate  the  effects  of 
uniform  compression  or  radial  translation,  if  indeed  such  would 
be  desirable,  but  ensured  that  the  effects  of  uniform  compression 
or  translation  did  not  dominate. 

It  may  be  noted  that  if  only  c„  ,  •»  and  i> ,  are 

adjusted,  then  the  deformation  is  one-dimensional,  i.e., there  is 
no  deformation  in  the  z  direction.  It  is  easily  verified  that 
for  such  a  deformation  all  of  the  finite  difference  equations 
immediately  reduce  to  the  well  known  second  order  one-dimensional 
difference  analog  for  «.  • I 

»  n 
a  •  - r  4.2 

i{  /}(*.-*')} 

where  we  note  that  p,  ■  pK  ,  p%  •  P+  etc.  It  is  therefore 
possible  to  check  the  error  of  the  one-dimensional  difference 
analog  directly. 

Results  are  shown  in  Fig.  2,  where  the  percentage  error 
is  plotted  vs.  the  coefficient  c„  .  Deformed  mesh  shapes  are 
also  shown  to  give  a  physical  picture  of  the  deformation.  It  is 
quite  clear  that  the  error  grows  to  very  large  proportions  as  the 
distortion  becomes  large. 
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Results  of  a  number  of  calculations  in  which  several 
of  the  quadratic  coefficients  were  varied  incrementally  are 
shown  in  Fig.  3  through  8  (for  <*  »  /  )  .  In  each  case  the  next 
increment  in  deformation  coefficients  led  to  inversion  of  one 
or  more  meshes,  i.e.,one  or  more  mesh  areas  became  negative. 
Errors  below  170  were  not  considered  significant.  Deformed  mesh 
shapes  are  again  shown  to  allow  the  deformation  to  be  visualized, 
lhe  errors  grow  very  rapidly  as  the  deformation  becomes  extreme. 
The  difference  in  error  of  the  various  finite  difference  equa¬ 
tions  is,  in  most  cases,  unimportant  compared  to  the  magnitude 
of  the  error.  However,  results  labeled  5-13  are  generally 
somewhat  better  than  these  labeled  1-4  (Taylor  Expansion  and 
Midpoint  Method).  One  exception  may  be  noted  in  Fig.  7  where 
method  10  (Eq.  3.22  using  Eq.  3.24)  is  much  better  than  all  the 
others.  In  all  cases  the  error  in  angle  did  not  exceed  a  few 
degrees  for  the  most  extreme  deformations,  and  was  quite  neg¬ 
ligible  over  most  of  the  range. 

Similar  results  for  ol  •  2.  are  shown  in  Fig.  9  through 
16,  the  errors  in  both  magnitude  and  angle  being  shown.  Defor¬ 
mations  correspond  to  those  in  Fig.  3,  4,  5  and  7  respectively. 
The  errors  in  magnitude  are  generally  smaller  than  the  corre¬ 
sponding  errors  for  <x  »  I .  None  of  the  errors  approach  zero 
as  the  deformation  is  reduced.  This  may  be  traced  to  the  ef¬ 
fect  of  radial  translation,  and  is  not  a  serious  concern.  It 
is  clear  from  the  results  that  no  one  finite  difference  equa¬ 
tion  is  clearly  superior,  and  errors  of  20-307.  in  magnitude  and 
20-30  degrees  in  angle  might  be  anticipated  under  certain  con¬ 
ditions  with  any  of  them. 

While  crossing  of  mass  points  could  be  observed  by 
suitable  adjustment  of  the  quadratic  coefficients,  this  occurred, 
for  the  combinations  used  above,  during  very  severe  distortions 
where  the  error  was  extremely  high.  The  effect  of  crossing  of 
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mass  points  is  somewhat  more  easily  illustrated  by  choosing  a 
different  deformed  mesh  shape  for  which  the  true  acceleration 
was  not  found.  Several  stages  in  crossing  of  mass  points  may 
be  distinguished. 

In  a  normal  mesh,  all  included 
interior  angles  are  less  than 
180°  as  shown. 

The  first  stage  of  crossing  oc¬ 
curs  when  one  included  interior 
angle  exceeds  180?  This  may  be 
termed  a  re-entrant  mesh. 

The  second  stage  of  crossing 
occurs  when  one  mesh  vertex 
crosses  an  opposite  side.  The 
mesh  then  takes  the  form  of  two 
triangles,  one  of  which  has  a 
negative  area.  This  may  be 
termed  a  schizoid  mesh. 

If  the  deformation  proceeds  far 
enough,  the  net  area  of  the  mesh 
becomes  negative.  This  may  be 
termed  an  inverted  mesh. 

The  effect  of  each  stage  was  investigated  by  solving 
the  finite  difference  equations  for  the  mesh  configurations 
shown  in  Fig.  17.  Configuration  (a)  corresponded  to  a  normal 
mesh.  The  acceleration  was  directed  to  the  original  position 
of  the  central  vertex,  for  «  *  I  ,  as  might  be  expected.  (For 
oc  •  2.  »  the  ret  radial  displacement  of  the  centroids  of  the 

meshes  leads  to  an  acceleration  directed  in  a  slightly  different 
direction) . 

In  configuration  (b)  the  upper  left-hand  mesh  becomes 


Normal 
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re-entrant.  All  finite-difference  equations  gave  the  expected 
acceleration  direction  except  10  (Eq.  3,22  using  Eq.  3.24)  which 
showed  an  error  near  180?  This  reversal  in  acceleration  would 
lead  to  an  accentuation  of  the  deformation,  and  hence  instabil¬ 
ity.  In  configuration  (c)  the  upper  left-hand  mesh  became 
schizoid.  Here  both  9  (Eq.  3.20  using  Eq.  3.21)  and  10  (Eq.  3.22 
using  Eq .  3.24)  showed  an  acceleration  reversal,  the  other 
finite-difference  equations  giving  the  correct  direction.  In 
configuration  (d)  the  upper  left-hand  mesh  has  become  inverted. 
All  equations  gave  an  acceleration  reversal.  (This  is  traceable 
in  part  to  the  fact  that  the  density  in  the  inverted  mesh  be¬ 
comes  negative.  The  computer  library  routine  for  raising  a 
number  to  a  floating  point  exponent  then  gave  a  negative  pres¬ 
sure,  from  equation  2.9,  which  was  sufficiently  large  to  out¬ 
weigh  the  contributions  of  the  negative  pressures  in  the  other 
expanded  meshes.) 

In  assessing  the  results  of  the  calculations  presented 
in  this  report,  it  seems  that  none  of  the  finite  difference 
equations  which  were  used  are  notably  superior  in  accuracy.  All 
except  9  and  10  (Eq.  3.20  using  Eq.  3.21  and  Eq.  3.22  using 
Eq.  3.24)  showed  no  reversals  in  direction  of  the  acceleration 
until  at  least  one  mesh  area  became  negative.  At  this  point 
all  methods  failed. 

While  the  present  calculations  are  hardly  extensive 
enough  to  warrant  drawing  general  conclusions,  it  does  appear 
that  any  of  the  methods  are  equally  satisfactory,  provided  a 
sufficiently  small  mesh  size  is  chosen.  It  seems  appropriate 
to  choose  the  method  which  requires  the  least  number  of  arith¬ 
metic  steps  for  solution.  Equations  3.20  seem  most  suitable 

for  this  reason. 
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Figure  3  Pure  Shear,  a  -  1 
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Figure  7  Shear  and  Rotation,  a  -  1 
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Figure  9  Pure  Shear,  a  •  2  Magnitude  Error 
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Figure  10  Pure  Shear,  a  ■  2  Angular  Error 
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Figure  13  Sheer,  Rotation  and  Compression,  a  -  2  Magnitude 
Error 
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